!
!***********************************************************************
!
      subroutine spline(x, y, n, yp1, ypn, y2, u) 
!-----------------------------------------------
!   M o d u l e s 
!-----------------------------------------------
      USE vast_kind_param, ONLY:  double 
      use modify_com_M 
 
!...Translated by Pacific-Sierra Research 77to90  4.3E  14:13:36   8/20/02  
!...Switches: -yf -x1             
      implicit none
!-----------------------------------------------
!   D u m m y   A r g u m e n t s
!-----------------------------------------------
      integer , intent(in) :: n 
      real(double) , intent(in) :: yp1 
      real(double) , intent(in) :: ypn 
      real(double) , intent(in) :: x(*) 
      real(double) , intent(in) :: y(*) 
      real(double) , intent(inout) :: y2(*) 
      real(double) , intent(inout) :: u(*) 
!-----------------------------------------------
!   L o c a l   V a r i a b l e s
!-----------------------------------------------
      integer :: i, k 
      real(double) :: sig, p, qn, un 
!-----------------------------------------------
!
!=======================================================================
!
      y2(1) = -0.5 
      u(1) = (3./(x(2)-x(1)))*((y(2)-y(1))/(x(2)-x(1))-yp1) 
      do i = 2, n - 1 
         sig = (x(i)-x(i-1))/(x(i+1)-x(i-1)) 
         p = sig*y2(i-1) + 2. 
         y2(i) = (sig - 1.)/p 
         u(i) = (6.*((y(i+1)-y(i))/(x(i+1)-x(i))-(y(i)-y(i-1))/(x(i)-x(i-1)))/(&
            x(i+1)-x(i-1))-sig*u(i-1))/p 
      end do 
      qn = 0.5 
      un = (3./(x(n)-x(n-1)))*(ypn - (y(n)-y(n-1))/(x(n)-x(n-1))) 
      y2(n) = (un - qn*u(n-1))/(qn*y2(n-1)+1.) 
      do k = n - 1, 1, -1 
         y2(k) = y2(k)*y2(k+1) + u(k) 
      end do 
      return  
      end subroutine spline 
